{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 6, demo 4\n",
    "\n",
    "Bayesian data analysis\n",
    "\n",
    "Posterior predictive checking  \n",
    "Light speed example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import preliz as pz\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "pz.style.use('preliz-doc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# data\n",
    "data_path = os.path.abspath(\n",
    "    os.path.join(\n",
    "        os.path.pardir,\n",
    "        'utilities_and_data',\n",
    "        'light.txt'\n",
    "    )\n",
    ")\n",
    "y = np.loadtxt(data_path)\n",
    "# sufficient statistics\n",
    "n = len(y)\n",
    "s2 = np.var(y, ddof=1)\n",
    "my = np.mean(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# tail area probabilities of marginal predictive distributions\n",
    "Ty = pz.StudentT(n-1, my, np.sqrt(s2*(1+1/n))).cdf(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA/MAAAGbCAYAAACIxMC9AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAR15JREFUeJzt3Xl8Dffi//F3ZCELsYQuoeJypalStLVrbdWWBrXTKqVFlaLVXnW73Kqq1r7WVpQqrZZWLdXa1yALqvY9SVXsJEgiOb8//M58M+ckkciJGF7PxyOPR2bmzJzPzPlMct7z+cxn3Gw2m00AAAAAAMAy8uV1AQAAAAAAQPYQ5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACLIcwDAAAAAGAxhHkAAAAAACyGMA8AAAAAgMUQ5gEAToKDg42fBg0auHTb48ePN21/4cKFLt0+bk2nTp1Mn0tMTExeF+mexWcBAMgKj7wuAAAg9wQHB5umV61apZIlS+ZRaVxj7969WrlypTFdrVo1Va9ePQ9LBAAAcPsR5gEATgICAozfixQpkoclcbZ3715NmDDBmO7duzdhHgAA3HMI8wAAJ5s2bcrrIgAAACAT3DMPAAAAAIDF0DIPAHCS9l77wMBArV692uk1K1eu1Ndff619+/bJ09NTlSpVUs+ePXX//ferYcOGxuuqVaumOXPmZPp+586d04QJE7R69WqdOXNGJUqU0PPPP6/evXvL29tbkrR161a98sorTutOmDDB1O3+xRdf1LBhw7K0n/v27dPs2bMVERGhU6dO6fr16ypUqJCKFSumChUqqEqVKgoNDZWPj48kKSYmxmnfZs2apTlz5uinn37S8ePH5ePjo5o1a6pv374KCgpK932vXbumn376SX/88Yf279+vy5cvy9fXVyEhIWrevLmaN2+ufPnSv95+/fp1LV26VMuWLdNff/2lCxcuyNvbW+XKldPzzz+v9u3by8vLK911o6OjNW7cOG3cuFFXrlxR6dKl1a5dO3Xs2DFLx+tmLl26pPnz52vNmjU6cuSIEhIS5O/vr4oVK6p169Zq1KiR6fXbt2/XK6+8otTUVElS2bJl9fPPPxvlT01NVceOHRUVFWWsM3HiRGM7W7Zs0fr167V79279888/On/+vK5evSpfX18FBQWpdu3a6tixo4oXL+5U1gYNGig2NtaY3r9/vxYtWqQ5c+bo8OHD8vf3V+PGjdWvXz/5+fkpMTFRU6ZM0eLFi/XPP/+oRIkSatKkiXr37q0CBQqYtj1w4EAtWrTImJ49e7aKFSum8ePHa9u2bUpISFBQUJDatWunDh06ZPhZZyYn9QAAcHcgzAMAsm3KlCkaNWqUad6GDRu0adMm9enTJ1vbOnz4sEaNGqXTp08b82JjYzV9+nTt27dP06dPl5ubm0vKndbGjRvVs2dPJScnm+afPXtWZ8+e1YEDB7Ro0SJVqVJF5cuXT3cbycnJeuONN7Ru3TpjXmJiopYtW6Z169Zp1qxZqlSpkmmdw4cPq1evXjp27Jhp/oULF7RlyxZt2bJFCxcu1MSJE1WoUCHTa+Li4tS7d2/t3LnTqRyRkZGKjIzUggULNG3aNN1///2m1+zdu1evvPKKLl26ZMzbv3+/Bg8erG3bthmB+lZFRkbqrbfeMn2OknTmzBmtWbNGa9asUZMmTfTFF18YIfPJJ5/Uq6++qq+//to4Nl999ZX69u0rSfr2229NQb5Vq1amCwKzZs3S2rVrncpy8eJF7dy5Uzt37tR3332nGTNmqEKFCpmW/9NPP9W3335rTF+7dk1z5sxRZGSkZs2apW7dumnXrl3G8tjYWE2bNk379+/XtGnTMt325s2bNXPmTCUmJhrz7Mc+IiJCI0eOzFYdz0k9AADcPehmDwDIlvDwcI0ePdppfoECBZSamqqxY8dma3vTp0/X6dOn5eHhIU9PT9OyjRs3av369ZIkT09PBQQEyM/Pz/QaHx8fBQQEGD+OyzMyatQoU5D38PBQ4cKFs9VKGhUVZQR5ew8Cu4SEBPXv31/Xrl0z5l24cEGvv/66U5B3LPO2bdv07rvvmuYlJSWpZ8+eTgHO19fXFAQPHDigN954Q0lJSaZ1+/XrZwryacv822+/KTIyMiu7nK4TJ06oR48epiDv5ubmtF/Lli1z6jXRr18/08USe0COiYkx1bOSJUtq0KBBGZbB09NThQsXlq+vr2n+hQsX9N5778lms2W6D/Yg7/g5/vXXX2rWrJkR5B1b4devX68NGzZkuu3JkycrMTFR+fPndwrtS5cu1fz58zNdP62c1AMAwN2FMA8AyJavvvrKFIwqVaqktWvXaseOHZo2bVqWw3Rar7/+usLDw7Vt2zan59rbg1LVqlW1adMm/fe//zUt79q1qzZt2mT8fPDBB1l6zwMHDhi/N2vWTOHh4dq6dat27dqllStXasiQIWrYsKHTBQZH1apV04YNGxQVFaVZs2aZWtNjYmK0ePFiY3rGjBmmrt316tXT+vXrFRERofXr16tq1arGsrVr15oGIly0aJH++usvY7pSpUpGCN+6daueeeYZY9mePXv0888/G9O//fab6QKCj4+Ppk6dqqioKG3evFk1atTIUcv82LFjTRcKWrVqpa1btyoiIkLLli1TmTJljGXz5s3TkSNHjGkvLy8NHz7cOM7Jycn64IMP9OGHH+rKlSuSpHz58umLL75wqlsvv/yyvv/+e0VFRWn37t3aunWrIiMjtWXLFjVt2tR43aFDh0yt6ukJDAzU0qVLtWPHDqcLKSdPnlSFChW0fv16RUVFqWXLlqblaXtmpMfT01PDhw9XZGSkwsPD9cILL5iWT5kyRSkpKZluwy4n9QAAcHchzAMAsuzq1asKCwszzfv000/1wAMPyM3NTU899ZQ6d+6crW2GhIRowIAB8vb2lo+Pj7p162ZaHh0dneNypydtC667u7sRZj09PVWqVCm1adNGkyZNMgVRR/ny5dPQoUNVokQJubm5qWbNmurSpYvpNWm7gS9dutT43cvLSyNGjNB9990nSbrvvvv03nvvmdZdsmRJuutK0rBhw4yy+fv766OPPspwXcew+dJLL+npp5+Wm5ubihUrps8++0zu7u4Z7mdmkpKStHLlSmO6RIkSGjJkiPz9/SXduA++d+/exvLU1FQtW7bMtI2HH35Yb731ljG9a9cubd682Zju1q2bnnjiCaf3rlu3rlJSUjR8+HC1b99ejRo1Up06dRQaGuq0z3v27Ml0P9566y2VK1dOkhQaGuq0fNCgQbrvvvuUL18+tW3b1rTsZnU0NDRUzZo1k4eHh/z8/DR48GAVLFjQWH7y5EnTxaXM5KQeAADuLtwzDwDIsuPHj+v69evGdEBAgB5++GHTa2rVqqWJEydmeZuOLfFFixY1TdtbZ12tfv36xiBlixYt0i+//KLAwECVLVtWwcHBqlmzpqpXr55pt/tSpUqpVKlSpnk1atTQuHHjjOnDhw9LutHtPiYmxpiflJSUbkBNa/fu3cbv+/fvNy1r0qRJlte1lyFtGdMqWbKkSpYsqePHj2e6zfQcO3bMdCtBXFycQkJCslw2u9dee01r165VRESEab5j0E9r8ODBmjt3bpbKeeHChUyX16xZ0/jdsQ56enrq8ccfN6YDAgJMy29WR9NuW7pxIalixYqmCxaHDh266XGTclYPAAB3F1rmAQBZFh8fb5p2DD0ZzcuM4wBdN+vW7iqDBg1S/fr1jenU1FRFR0dr7dq1mjJlirp06aJmzZqZArij9Pa1SJEipumEhARJzscuK86fP2/8nt31ExISjPul7WXIqIwZzcuKy5cvZ3ud9IJ1vnz51Lp1a6f5oaGh6Y7Kvnr16iwHeUmmi1DpKVasmPG7Yx0sUqSI6X707A7ImJ16cjM5qQcAgLsLLfMAgCxzvGc5vVB27ty5bG3Tw8P8ryg3Rq5PT6FChTR58mQdO3ZMmzZt0sGDB3XixAn9+eefxv3fBw8e1NChQzVp0qR0t5E2bGc0z96d3/HY2Qfcy4y9q7p9ffvxtnePvxn7fdiOg8JlpdxZlba7uHTj9gHHUfgdOZZHulFvHJ+QIEmTJk3Sc889p5IlS5rm//7776bp5557Tm+//bYCAwPl4eGh+fPn6+OPP87qbjjVw7RyeoEpK8c7q2NN5KQeAADuLoR5AECWBQUFydPT0xgF/vTp0zp27Jjpeeppuw7nBsdu7zkNKkFBQabyX7t2TS+88IJxH/TWrVszXDc6OlqxsbEKDAw05jmOKVC2bFlJNwJsyZIljZZ+Dw8PrVixItMQl3ZQuuDgYKMsNptN3333nUqXLp3puvZjVbZsWe3du9dUxjp16hjTMTExmfZAyExQUJAKFChgdLUvXry4Vq5cmentCekNtve///3PaTR8m82mhIQEDRw4ULNnzzZtMy4uzrR+r169TMcjJ6Pzu1pYWJjpPvyEhAT9+eefptfY68nN5KQeAADuLvx1BwBkWYECBUz3W9tsNn344YeKi4uTzWbThg0b9M033+RqGRzD786dO2+pG/Hbb7+tqVOnas+ePab1jxw5YhqZPbPu2SkpKRo0aJBOnz4tm82mLVu2aNasWabXpO3K//zzzxu/X7t2TX369DHdA52SkqIjR45o3rx56tKli2kk/LTrSlLfvn21Y8cO48kCNptN0dHRWrRokd544w1NmTLFeG29evVM686dO1fr1q2TzWbT2bNn9cEHH9zyRREvLy/TuAexsbEaMGCAaVC4pKQk7du3T7NmzVKbNm0UHh5u2sbPP/+sFStWGNMtWrTQyy+/bExv375dM2fONK3jWA9+/fVXpaSkKDk5Wd9++63p2OW1X375RUuWLFFKSori4+P10UcfmW5PePDBBxUcHJylbeWkHgAA7i60zAPAPaR169YZjlretWtXp5Hk09OzZ09t3LjRCA/btm1T3bp15e3tratXr+Z6N3nH0LN582Y9/vjjRtfukSNHOg3wlp7Dhw9r6dKlGjlypPLly6eCBQvq+vXrTvcuV6pUKcNt5MuXz2jltu9/WiVLljS1yHbr1k1Lly7V33//bZS9WbNm8vLyko+Pj+Lj400XD5o1a2b83qpVK/3www/GqOx79+5Vu3btjBHSExISjB4TkvTII48Yvz/77LOaMGGC8Xi6K1euqHv37umW+Vb069dPGzZsMALq0qVLtXTpUnl7e8vLy0vx8fGmiwVpH2148uRJDRkyxJguXry4Bg0aJC8vL61fv94YlG/MmDGqW7eu8Uz6unXrmi4ATJs2TbNnz1ZqaqqSk5NNvQXyWkpKit555x0NGjRIycnJTj0TunfvnuXW85zUAwDA3YWWeQC4h5w/f15nzpxJ9yero8Y/8cQT6t+/v9P8q1evysPDQwMGDDDNd3W4L1WqlJ566inTvKSkJGM/bqWVPjU1VRcvXnQK8v7+/ho0aFCG61WtWlXPPfecJDmFYl9fX40ePVoFChQw5hUpUkTTp0/Xv/71L6fyX7hwwakXgI+Pj/G7l5eXpkyZoipVqphec/36dV24cMEU4NJbd/To0U73stvLXLt2baftZkfp0qU1depUlShRwmn7Fy9eNAV5d3d345jYbDYNHDjQ1Er9v//9T/7+/vL29tbQoUON+pOUlKR3333X+HybN2+uypUrm94vMTFRycnJCgwM1JtvvnnL++Nq/fv3l4+PjxITE52CfNOmTdW+ffssbysn9QAAcHehZR4AkG09evRQ2bJlNX36dO3bt0+enp6qUqWKevXq5RRI0w7i5iqjR4/WuHHjtHr1av3zzz9OASYrhg4dqo0bN2r79u2Kjo7WuXPnlJCQIG9vb5UuXVq1atVSp06djOfApydfvnwaPXq0nnzySf3www86duyYvL29VatWLfXt29d0L75d2bJl9fPPP2vx4sX6/ffftXfvXl24cEH58uVTsWLFVK5cOVWvXl2NGzfWQw89ZFq3RIkSmjt3rlasWKFly5bpzz//1Llz52Sz2VSkSBGVLVtWTzzxhBo1auT0yMBHHnlECxcu1NixY7Vp0ybFx8froYceUvPmzdW1a1e9+uqr2T6GaVWtWlXLly/Xjz/+qNWrV+vgwYO6dOmSPDw8VLx4ceNxf88884xxTL/55hvTGANNmjRRo0aNjOknnnhCL7/8subMmSNJ2rdvn8aPH6933nlHXl5emjlzpsaPH6/ly5frzJkzKlasmOrXr6++fftqzZo1OdofV3rssce0cOFCjRs3TmFhYYqPj1dQUJDat2+vDh06ZPuCV07qAQDg7uFmS9vXDQCAHJoyZYppVPKePXum25JvRTExMWrYsKExXa1aNSNoAnYDBw7UokWLjOnZs2erevXqeVgiAMDdiG72AIBsGz16tOm+ebuNGzc6DbjVuHHj21k0AACAewLd7AEA2RYZGanJkyerSJEiCg4Olru7u6Kjo3XixAnT61q0aKEKFSrkUSkBAADuXoR5AMAtO3/+vNNz1e1CQ0P16aef3uYSAQAA3BsI8wCAbOvRo4fKlCmjqKgonT59WpcvX1b+/Pn1wAMPqHLlymrZsqUef/zxvC4mAADAXYsB8AAAAAAAsBgGwAMAAAAAwGII8wAAAAAAWAxhHgAAAAAAi2EAPOAu1qBBA8XGxhrT+/fvN37funWrXnnlFWP6xRdf1LBhw25r+bLqXtiPu11CQoKmTp2qP/74QzExMUpMTDSWbd++XYUKFcrD0lnDwoUL9f777xvTvXv3Vp8+ffKwRP/nXq7bd4I7uW7cKzgHAOQFwjwAl4mJidGiRYuM6ZCQEDVq1CgPS5T77sV9zi6bzabXXntNkZGReV0UwCVmzZqly5cvG9MEZwBAXiDMA/coT09PBQQEGNN+fn453mZsbKwmTJhgTL/44osuCbZFihQxteTeSbKzz3fyfuSmLVu2OAV5X19feXt7S5Ly5eOOr6woUKCA6Zz18fHJw9Lc22bPnm1qhc3rME/dAIB7E2EeuEdVrVpVmzZtyutiZMlPP/2U10VwibtlP7Jrz549pumOHTvq448/zqPSWFeTJk3UpEmTvC4G7kDUDQC4N9EcAgDIVdeuXTNNV6xYMY9KAgAAcPegZR6wuH379mncuHEKDw9XcnKy/v3vf6tz585q2rRppuvdbOC41NRULV68WL/++qv279+vCxcuyMPDQ4ULF9aDDz6oypUrq3bt2qpdu7ZiYmLUsGFDp/dYtGiR6X7yatWqac6cOZLSHyxo+fLlmj17tg4cOKD4+HitWrVKJUuWzPbAQsnJyZoxY4Z+/vlnxcTEyN/fX/Xr19dbb72l4sWLm147fvx4Uzf5zz//XC1btjSmHffNvg+u2mdHV65c0Q8//KBVq1bp4MGDunz5snx9fRUUFKSnnnpKHTt2VNGiRZ3WS2/bGzZs0IwZM/Tnn38qOTlZ5cuXV/fu3fXMM89kevwykt2yOQ7KZff+++8b89Men8yk9zlVqlRJY8eO1datW5WSkqLKlSurb9++qlSpkiRp7dq1mjp1qvbu3SsPDw89/vjjevvtt1W+fHnTtpOSkrRgwQLt3r1b+/bt09mzZ3XhwgVJUuHChRUSEqImTZooNDTU6ZaA9M6jDz/8UF999ZVWrFihkydPqkqVKqZ9DA8P16RJk7Rz507ZbDaFhISoa9euatiwoYKDg43XBQYGavXq1RkeT8dBzjp16qRt27YZ06tWrdK5c+c0efJkRURE6MqVKypTpoxeeukltWvXzukYb9myRevXr9fu3bv1zz//6Pz587p69arxGdeuXVsdO3Z0OodyKr3PtlatWho3bpzWr1+vixcv6sEHH1SLFi3UrVs3eXl5pbudw4cPa+7cudq2bZv+/vtvJSUlqUiRIqpYsaJCQ0P17LPPpntLx759+zR79mxFRETo1KlTun79ugoVKqRixYqpQoUKqlKlikJDQ+Xj4+N0jO3Sfm6SjL9ddnFxcZo7d642bNig6OhoXb16VUWLFlWVKlX08ssv68knn3TaZnqfd8eOHTVu3DitW7dOcXFxatasmYYNG5blAfB27dqlefPmKTIyUnFxcUpJSVFAQIAee+wxtW7dWrVr185ROW5m4MCBpr+Ps2fPVrFixTR+/Hht27ZNCQkJCgoKUrt27dShQ4ds34LTtWtXU4+zP/74Qw899JDpNUeOHNHzzz9vTNepU0dff/21pNw7B7Ly+WR27ttdunRJ8+fP15o1a3TkyBElJCTI399fFStWVOvWrTO8xSs7dRyAtRDmAQvbsGGDevXqpaSkJGPezp079fbbb2v37t052vagQYNMX7qkGyH56tWrOnnypCIiIhQWFpbul79bMXr0aE2ePDnH20lMTNSrr76q7du3G/NOnz6tH374QWvXrtV3332nUqVK5fh9csO+ffvUq1cvUyiXpIsXL2rnzp3auXOn5syZozFjxqhmzZqZbmvs2LGaNGmSad6uXbvUu3dvDR8+XM2aNcuzsrlCRESEPvnkE1Or/8aNG7V9+3bNnDlTO3bs0JdffmlaZ82aNdq+fbt+/PFHlSlTxpgfHx+vwYMHp/s+p06d0qlTp7R27VotXrxYkydPlqenZ4blunTpktq3b68DBw6ku3zx4sX6z3/+o9TUVGNeeHi4wsPD9fbbb2dp37NqwYIFmjp1qum99u/fr48++kjnz59Xz549Ta+fNWuW1q5d67SdtJ/xd999pxkzZqhChQouLWta+/fv1xdffGFcUJGkY8eOacyYMdq8ebOmT5+u/Pnzm9aZOXOmhg8frpSUFNP8uLg4rVq1SqtWrVKNGjU0btw4+fv7G8s3btyonj17Kjk52bTe2bNndfbsWR04cECLFi1SlSpVnC4CZdXKlSv13nvvKSEhwTT/1KlT+u233/Tbb7+pS5cuGjhwoNzc3DLczt9//60WLVooLi4u22Ww2Wz64osvNHPmTKdlsbGxio2N1bJly9S0aVMNGzYswwsmOS2Ho82bN2vmzJmmsUT279+vwYMHKyIiQiNHjsz0mDhq3ry5KcwvX75cPXr0ML1m+fLlpukWLVoYv98p50B6IiMj9dZbb+n06dOm+WfOnNGaNWu0Zs0aNWnSRF988YXp87sddRxA3qGbPWBR586d04ABA0xBXpIxqNiMGTN08uTJW9q2/Z+743Z9fX3Tfb27u7sCAgJMX5IlKX/+/AoICDB+HJenZQ/y+fPnz1HrwG+//abt27fLzc3N6Qt/XFycBgwYIJvNdsvbt3PFPqd17tw5vf76605h2f552l24cEFvvvmmjh49mun27EG+QIECTsvSCz25UTb7oFyOn6efn1+2j4+jH3/8UdeuXVP+/PlNX/YTExPVv39/jRgxwihDWvHx8Ro/fnym2/bx8VGRIkWcAs3GjRvTDUNprVq1ygjyhQoVkru7u7Hs+PHj+vDDD03hWvq/4zhq1KhMt51dkydPVmpqqtN5IN2oHxcvXsxwXU9PTxUuXNjpnL9w4YLee+89l5xDGZk1a5YuXLggT09PeXiY2xy2bdvm9PktXrxYw4YNM9VpNzc3p88+LCzM6YLJqFGjTCHH3vsooxZhf39/BQQEOC1Pe84HBAQYn3tUVJT69etnCvL58uVzOq6zZs3SjBkz0n1Pu4ULFyouLk5ubm4qVKhQtkLu5MmTnequu7u7Ux1funSphgwZkmvlSK9ciYmJTuexvSzz58/P1vYaN25s+nvjGNylG/8j7Hx9fTPsqZSX54CjEydOqEePHqYg7+bm5jRw7bJly5x6SGS3jgOwFs5kwKK+//57U8tV8eLF9f3332vHjh1auXKlgoODnUJDVjm2Kn711VfasWOHIiMjFRERoQULFqh37956+OGHJUkPPPCANm3a5PQlu0mTJtq0aZPxk7YrrSMvLy99+eWXioyMVFRUlBYvXqwiRYpku+ypqal6/vnntW3bNkVFRWnUqFGmltQdO3Zo8+bN2d6uI1fsc1rTp083tXSVLl1av/zyi3bs2KHVq1frscceM5YlJCRo7NixmW7P399fs2bN0o4dO/Tzzz+rcOHCxrK4uDjt27cvS+XKSdnsx6Jr166m7f33v//N9vFx5Obmpk8++USRkZH6/fffTV+4T506JQ8PD02cOFE7duzQN998Y1p3/fr1pi/iPj4+GjlypP744w/t2bNHUVFRCgsL065du7Rw4ULTKOGOF7nS88gjj2jZsmXavn27duzYof/85z+SpK+//trUk6B06dJasmSJduzYoQULFqhEiRK3dCwy4uXlpdGjRysqKkqrVq0ydTdOTEzU1q1bTa9/+eWX9f333ysqKkq7d+/W1q1bFRkZqS1btphu2zl06JB27drl0rI6GjhwoPH3plu3bqZl3377rfFYuKSkJOPCjV2bNm0UHh6uyMhITZw40RTuNm7cqHXr1hnTaf/WNWvWTOHh4dq6dat27dqllStXasiQIWrYsKHxN2TChAnatGmTHnjgAdN7pj3n0y7/4osvTEGqZ8+eioyMVGRkpBYsWKBixYoZyyZMmJDpBRZJql27ttHDJDIyUl26dMn09dKNi3FTpkwxzevdu7dRjk8//dR00emHH37QwYMHXV6O9Hh6emr48OGKjIxUeHi4XnjhBdPyKVOmZOvCo7e3t5599lljeu/evaYLn4cPHzZ95s8++6zpos+ddA6kNXbsWF26dMmYbtWqlbZu3aqIiAgtW7bM1NNo3rx5OnLkiDGd3ToOwFoI84BFpf1CKklvvvmmKleuLEkqVaqUPvroo1vetmNLar58+Yzw4+fnp0qVKqlPnz4aOnToLb+Ho1deeUXNmzc3WuKCg4Mz7AmQGT8/Pw0ZMsRoFW3atKlTl/L0ulHmtbStRdKN2xzsF0sCAwOdWsvWrFnj1CsjrTfffFM1a9aUm5ubQkJCnO7vj46OzrOyuUL16tXVvn17eXh46KGHHjLqvl1oaKgaNWokNzc31ahRQ0FBQcayy5cv6/z588Z0gQIF1LhxY23dulX9+/dXaGio6tWrpzp16qh79+6mgHX06FGnAf3ScnNz05dffqmyZctKuhGoH330UUk3jktaAwcO1L///W9JUqVKldSvX79bORQZ6tChg5o0aSJ3d3eVLFnS1J1Ycq4DdevWVUpKioYPH6727durUaNGqlOnjkJDQ53+3jg+ocCVqlevrldffVVeXl4qUKCA3n33XZUuXdpYfvXqVeNCRFRUlE6dOmUsK1GihD7++GP5+fnJ3d1djRo1Uvv27U3bd2yZtXN3dzcugHp6eqpUqVJq06aNJk2aZApLWXXy5ElFRUUZ04899pj69+9v9MSoVKmSOnfubCy/cuWKUx1Jy9vbWyNGjDAuFPj4+BjnYWbWrl2rq1evGtOPPvqo+vTpowIFCsjT01Nt27Y13Wtts9m0YsUKl5cjPaGhoWrWrJk8PDzk5+enwYMHq2DBgsbykydPZnjLSkaaN29umk7bOp9ZF3vpzjkH0kpKStLKlSuN6RIlSmjIkCFGr6ayZcuqd+/exvLU1FQtW7bMmM7NOg4g73HPPGBRhw8fNk3XqFHDNP3444/L09PT6T65rHjyySdVsGBBo/WrR48e8vb2VpkyZfSvf/1LFStWVIMGDZwGFsoJxy9gt6pSpUpOXQ9r1Khheiyc47HLawkJCU5d2B3vOy9fvrwCAgJ05swZSTdGiD9+/LgRBh01aNDANJ22BVCS6cv97S6bKzjWd8deHNWrVzdNBwQE6NixY8Z02v0/deqUOnfufNNbF6QbQefixYvp3r4g3WiVT2+/L1++7HSPseN4E7Vq1brp+2dHduvA4MGDNXfu3CxtO22vIFdzrF9ubm6qVq2ajh8/bsw7dOiQGjVq5NSC/OSTTzq1MNasWdPUhT1tOKxfv77R22LRokX65ZdfFBgYqLJlyyo4OFg1a9ZU9erVb6lLsmPvl507dzoNlOdo9+7dTgHT7qmnnkp38MubcTxG6dWzmjVrmgJ8ZgH6VsuRHsfP2tfXVxUrVjT1njp06JBCQkK0bNkyffbZZ+lup2vXrkYPjho1auiBBx4wbjNbvny5evXqZfxuFxgYqGrVqpm2c6ecA2kdO3bMdAExLi5OISEhma6Tdsyc3KzjAPIeZy5gUY6DKTmGGTc3N1PX6uwoWLCgJkyYYOpKevXqVe3Zs0dLlizR559/rsaNG2vw4MEuu28wMDDQJdtJ70um47FxPHZpOe7P9evXXVKuzMTHx5umfX19073P2XHf7Bdb0nP//febph0DTlY/t9womyuk7fouOe+f43LHL6pp93/IkCFZCvJ2mV0gSzuCeVqOdS694+gYtnMqO3Vg9erVWQ4xUu6eF9k5hx3rWXrrOs5LW6cHDRqk+vXrG9OpqamKjo7W2rVrNWXKFHXp0kXNmjVTTExMtvfjVs6BtD1GHN3q38icHiNXlSM92fmsr127pjNnzqT7c+XKFeP1bm5upt5YBw4c0OHDh3Xo0CEdOnTImN+sWTPTffq38xxw/Pub2d+UW6lHaS805GYdB5D3aJkHLMrX19d0D9358+dN4d1ms+Wo5aBGjRpauXKltm/frqioKB09elRHjhzRnj17lJqaKpvNprlz56pWrVoZPg4nO26lS3160vsy7Dgvbcu946BLjl+q/vnnH5eUKzOOPQkSEhKMQaHSOnfunGk6bXdUR47B7VYHqcqNsrlC2nt80+M4cFpGkpKSTF2bPTw89PHHH6tx48bG+dSuXTvt2LEjS9vLaPBGx+N45coVJSUlmQYgczyGOeV4DDKrA7///rtp+rnnntPbb7+twMBAeXh4aP78+fr4449dWr6MZOccdqxn6R1Dx3lpP4tChQpp8uTJOnbsmDZt2qSDBw/qxIkT+vPPP42/rwcPHtTQoUOdng5xM45ly2wQ0bSvycitDgya02PkqnKkJ7t/r7OqefPmpnECli9f7jSGjGMPiNw8B3Lyf8bx8/Py8lKhQoUyfb+09Sw36ziAvEeYByyqbNmypvsxw8LCTPe8RURE3FIX+7Q8PDxUs2ZNU1fIjRs3mgakCgsLM8K8Y+tndgYucpVdu3YpISHB9GUmLCzM9Jp//etfxu+O3aUdH/uT2T2skmv22dfXV4GBgabu7Fu2bFG9evWM6QMHDhjd2KUb5U57H3FuuZPL5grnz583nSfBwcFq27atMX358uWbDgaWFX5+frrvvvuM+7ttNpvCw8NNXZ7TPlLrdnO8BaBXr16mzzAyMvK2lSUsLMz02Dybzeb0bHf7mASOtzRs375dycnJpotZW7ZsMb0mvcdvBQUFmcZVuHbtml544QVjXAHHwQIdw1lKSorTBSbHLvWPPPKIvvvuO6f3Tis3Rkh3PEabN2/WO++8Y5qXlWOUG8LCwhQaGmpMJyQk6M8//zS9xv5Zt2zZUi1btszSdsuWLatHH33U6G6+fPly07GtXLmy6fOWcvccyMn/maCgIBUoUMDoal+8eHGtXLky027x6Q1+m906DsAa6GYPWFTaMCVJEydO1M6dOyXdGNgqo+dmZ8WBAwfUq1cv/frrr6bH2yUnJzs9vz5tV0PHFoQ9e/Zk2qU9N1y+fFkfffSRLl++rJSUFC1btkyLFy82vSZtl0PHZ84vWrRIsbGxSklJ0fLlyzVv3rxM389V+5x2BGZJ+vzzz417bmNjY/XBBx+YlterVy/TZ0G70p1ctpzy8/MzBbMjR44Yo1THxcWpf//+LqvDaeudJA0dOtS4j3/Xrl0aM2aMS97nVji2fv76669KSUlRcnKyvv32W6dzKDdt2bJFs2bNUlJSkq5du6bhw4eb7pf39vY2xkSoUqWK6SkAcXFx+uSTTxQfH6/U1FStXLnS6fFmzz33nPH722+/ralTp2rPnj2mQRuPHDli6vnk2KXa8bx3vNggSQ8++KDpSQ8RERH67LPPTKHx2rVr2rVrl7766is1bdrUaXwKV6hXr54pTO7evVvjx4/XtWvXlJycrB9++ME0wJqbm5saN27s8nKk55dfftGSJUuUkpKi+Ph442+33YMPPnjTcQYykrbl/dChQ6axUtIboyU3zwHH/zOrVq3S3r17ZbPZFBYWpokTJ2a4rpeXl2n8i9jYWA0YMMA0gGVSUpL27dunWbNmGU9zsLvVOg7AGmiZByyqbdu2mjlzptGV/vTp02rbtq28vb2zPLhZRlJSUrRq1SqtWrVK0o0vE76+voqPj3dq7U/7ZTUoKEheXl7GF4ZDhw6pWrVq8vf3l5ubm959990MB3dypSVLlmjp0qXKnz+/08jjlStXNvU0qFmzpvLnz6/ExERJN74oNWjQwNQSkhlX7fPrr7+uX3/91WixOXbsmJo3by4fHx/T/aDSjW6uffv2vek2XeVOLltO+fr6qkqVKkar29WrV9WmTRv5+fkZ9w1ntS7cTNeuXfXzzz8b2zp48KCeffZZ45zNyfO6c6pu3bqmAdCmTZum2bNnKzU1VcnJyS47BlmRL18+ff755xoxYoRsNptTyHj55ZeNMO3l5aUBAwbovffeM5YvWLBAP/74Y7rnf506dfT0008b04cPH9bSpUs1cuRI5cuXTwULFtT169edLuBUqlTJNB0cHKy9e/ca0126dJG/v788PT1VoUIFTZ06VdKNJxa88sorxt/N2bNna/bs2fLx8ZGHh4cuX76c688rL1q0qHr27Gm6WDRhwgR99dVXcnd3d3ryRJs2bW5by3xKSoreeecdDRo0SMnJyU4tyt27d7/lgdmaNm3q9GhA6cYtSE2aNHF6fW6eAyEhIaaeOZcuXVKLFi2y/P+6X79+2rBhg3GhY+nSpVq6dKm8vb3l5eWl+Ph4U6+wtHXqVus4AGugZR6wqKJFi2rEiBFO90bbvxi8+OKLevDBB13yXklJSU7dkaUbo3Gn7SJZoEABtWrVyvSa69ev6+zZszpz5kyOLzJkRcOGDRUSEiKbzeb0xatEiRIaPny4KTQVKlTI9Fgfu2vXrsnNzU2vvfZapu/nqn0uWrSopk+f7jS4lGNYLly4sCZNmmS6VSC33cllc4X333/fqRusPch36NDBZV9yS5curcGDBzuFE3sdsT+P3u52hvvmzZs7Pd4vMTFRycnJCgwM1JtvvnnbytK9e3cFBAQoOTnZKchXq1ZNffr0Mc1r3ry5Bg4caOrmnt75X6NGDY0aNSrD901NTdXFixedQo6/v78GDRpkmte+fXunz/HixYs6c+aMaaySqlWrasyYMU4t+VeuXNGlS5dMocvT0zPXnvXds2dPp2fBp6SkOAX5pk2b6sMPP8yVMqSnf//+8vHxUWJiolOQb9q0qdNjBbOjaNGiqlu3rtP8+vXrpzs4bG6eA+7u7nrvvfeczmn7uf/6669nun7p0qU1depUUy8U+/oXL140BXl3d/cMn7aRnToOwBpomQcsrG7duvrxxx81duxYhYeHKykpSWXLllW7du3Utm1bp2eLZ1XZsmU1YcIEbd26VTt37tTp06d17tw5paamqnDhwgoODlaTJk1Mz4W3++9//6uAgAAtXbpUMTExuf68cUeFChXSyJEjNXnyZC1fvlwnT56Uv7+/6tevr7feekvFixd3Wqd79+4qUqSI5syZo6NHj8rb21tVqlRRz549Vbx4cU2fPj3T93TVPj/88MP69ddftWDBAq1atUoHDhxQfHy88VjAp59+Wh07dnTZY6HulrLlVKVKlTRv3jzjPEpJSVGZMmXUoUMHtW3bVp06dXLZezVv3lyBgYGm22IqVKig1157TcHBwRo2bJjxWvtzpG8HLy8vzZw5U+PHj9fy5ct15swZFStWTPXr11ffvn1vOnaEK5UuXVqLFi3SuHHjtG7dOp0/f16BgYFq3ry5XnvttXRv4Xj11Vf11FNPae7cuQoLC9PJkyeVnJyswoUL69FHH1WzZs303HPPOQXwoUOHauPGjdq+fbuio6N17tw5JSQkyNvbW6VLl1atWrXUqVMn3Xfffab1qlSpomnTphndl+Pj4zNsYW/UqJFWrFih+fPna8OGDTp69KguX76s/Pnz67777lNISIhq1aqlZ5555pafQHIzbm5uev/999W0aVPNmzdP4eHhiouLU2pqqooVK6bKlSurVatW6Ybf3PTYY49p4cKFGjdunMLCwhQfH6+goCC1b99eHTp0yPEFrRYtWmj16tWmeRk9BjW3z4EXXnhBBQoU0NSpU7V//365u7vr0UcfVdeuXVWvXj1NmzYt0/WrVq2q5cuX68cff9Tq1at18OBBXbp0SR4eHipevLjxmLlnnnnGVF9vtY4DsAY3W2737wIAADe1ZMkS08BkoaGhGjFiRB6W6PYYP368JkyYYEx//vnnWR7oDNYycOBA45nn0o3bDuzjHwAAso9u9gAA3CZff/21VqxY4dR9fNeuXU7B3XHgQQAAgLToZg8AwG1y8OBBffnll/Lz81NISIjy58+vv//+W0ePHjV1065WrZrxyEcAAID0EOYBALjN4uPjtX379nSX1alTR6NGjcrT0e0BAMCdjzAPAMBt0qFDB/n5+SkiIkKnTp3SpUuX5OnpqRIlSqhixYoKDQ3VU089RZAHAAA3xQB4AAAAAABYzG0ZAM9ms+nKlSsZPrYFAAAAAABk3W0J81evXlXnzp119erV2/F2AAAAAADc1Xg0HQAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACLIcwDAAAAAGAxhHkAAAAAACyGMA8AAAAAgMUQ5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACLIcwDAAAAAGAxhHkAAAAAACzGI68LAADAncxtxPi8LkKO2Ab0yesiAACAXEDLPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACLIcwDAAAAAGAxhHkAAAAAACyGMA8AAAAAgMUQ5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACLIcwDAAAAAGAxhHkAAAAAACyGMA8AAAAAgMUQ5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACLIcwDAAAAAGAxhHkAAAAAACyGMA8AAAAAgMUQ5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACLIcwDAAAAAGAxhHkAAAAAACyGMA8AAAAAgMUQ5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACLIcwDAAAAAGAxhHkAAAAAACyGMA8AAAAAgMUQ5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACLIcwDAAAAAGAxhHkAAAAAACyGMA8AAAAAgMUQ5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACLIcwDAAAAAGAxhHkAAAAAACyGMA8AAAAAgMUQ5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYjzyugAAAAB3K7cR4/O6CDliG9Anr4sAAMgALfMAAAAAAFgMYR4AAAAAAIshzAMAAAAAYDGEeQAAAAAALIYwDwAAAACAxRDmAQAAAACwGMI8AAAAAAAWQ5gHAAAAAMBiCPMAAAAAAFgMYR4AAAAAAIshzAMAAAAAYDGEeQAAAAAALIYwDwAAAACAxRDmAQAAAACwGMI8AAAAAAAWQ5gHAAAAAMBiCPMAAAAAAFgMYR4AAAAAAIshzAMAAAAAYDGEeQAAAAAALIYwDwAAAACAxRDmAQAAAACwGMI8AAAAAAAWQ5gHAAAAAMBiCPMAAAAAAFgMYR4AAAAAAIshzAMAAAAAYDGEeQAAAAAALIYwDwAAAACAxXjkdQEAAAAy4jZifF4XARZm9fpjG9Anr4sA4A5GyzwAAAAAABZDmAcAAAAAwGII8wAAAAAAWAxhHgAAAAAAiyHMAwAAAABgMYR5AAAAAAAshjAPAAAAAIDFEOYBAAAAALAYwjwAAAAAABZDmAcAAAAAwGII8wAAAAAAWAxhHgAAAAAAiyHMAwAAAABgMYR5AAAAAAAshjAPAAAAAIDFEOYBAAAAALAYwjwAAAAAABZDmAcAAAAAwGII8wAAAAAAWAxhHgAAAAAAiyHMAwAAAABgMYR5AAAAAAAshjAPAAAAAIDFEOYBAAAAALAYwjwAAAAAABZDmAcAAAAAwGII8wAAAAAAWAxhHgAAAAAAi/HI6wIAAIDc4zZifF4XAQAA5AJa5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACLIcwDAAAAAGAxhHkAAAAAACyGMA8AAAAAgMUQ5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACLIcwDAAAAAGAxhHkAAAAAACyGMA8AAAAAgMUQ5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACLIcwDAAAAAGAxhHkAAAAAACyGMA8AAAAAgMUQ5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACL8cjrAgAA7m5uI8bndREAAHnA6n//bQP65HURgEzRMg8AAAAAgMUQ5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACLIcwDAAAAAGAxhHkAAAAAACyGMA8AAAAAgMUQ5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACLIcwDAAAAAGAxhHkAAAAAACyGMA8AAAAAgMUQ5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACLIcwDAAAAAGAxhHkAAAAAACyGMA8AAAAAgMUQ5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYjzyugAAkNvcRozP6yLkiG1An7wuAgAAAO4wtMwDAAAAAGAxhHkAAAAAACyGMA8AAAAAgMUQ5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACLIcwDAAAAAGAxhHkAAAAAACyGMA8AAAAAgMUQ5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACLIcwDAAAAAGAxhHkAAAAAACyGMA8AAAAAgMUQ5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACLIcwDAAAAAGAxhHkAAAAAACyGMA8AAAAAgMUQ5gEAAAAAsBjCPAAAAAAAFuOR1wUAAGTObcT4vC4CAAAA7jC0zAMAAAAAYDGEeQAAAAAALIYwDwAAAACAxRDmAQAAAACwGMI8AAAAAAAWQ5gHAAAAAMBiCPMAAAAAAFgMYR4AAAAAAIshzAMAAAAAYDGEeQAAAAAALIYwDwAAAACAxRDmAQAAAACwGMI8AAAAAAAWQ5gHAAAAAMBiCPMAAAAAAFgMYR4AAAAAAIshzAMAAAAAYDGEeQAAAAAALIYwDwAAAACAxRDmAQAAAACwGMI8AAAAAAAWQ5gHAAAAAMBiCPMAAAAAAFgMYR4AAAAAAIshzAMAAAAAYDGEeQAAAAAALIYwDwAAAACAxRDmAQAAAACwGI/b+Wb+46bousdtfUuXsQ3ok9dFuKe5jRif10XIEavXH6sffwDAreHvP+5lVq//Vv/+iZujZR4AAAAAAIshzAMAAAAAYDGEeQAAAAAALIYwDwAAAACAxRDmAQAAAACwGMI8AAAAAAAWQ5gHAAAAAMBiCPMAAAAAAFgMYR4AAAAAAIshzAMAAAAAYDGEeQAAAAAALIYwDwAAAACAxRDmAQAAAACwGMI8AAAAAAAWQ5gHAAAAAMBiCPMAAAAAAFgMYR4AAAAAAIshzAMAAAAAYDGEeQAAAAAALIYwDwAAAACAxRDmAQAAAACwGMI8AAAAAAAWQ5gHAAAAAMBiCPMAAAAAAFgMYR4AAAAAAIshzAMAAAAAYDGEeQAAAAAALIYwDwAAAACAxRDmAQAAAACwGI+8LgAAAAAAZ24jxud1EQDcwWiZBwAAAADAYgjzAAAAAABYDGEeAAAAAACLIcwDAAAAAGAxhHkAAAAAACyGMA8AAAAAgMUQ5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACLIcwDAAAAAGAxhHkAAAAAACyGMA8AAAAAgMUQ5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACLIcwDAAAAAGAxhHkAAAAAACyGMA8AAAAAgMUQ5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACL8bgdb2Kz2W68Wcr12/F2ueLKlSt5XYR7msd169Ydyfr1x+rHHwAA4F5j9e+fdwNvb2+5ubnl2vbdbPaknYtOnjypt956K7ffBgAAAACAO8K0adNUuHDhXNv+bWmZ9/X11Z49e7RmzRr5+fndjrcE7hgJCQl66qmntH79evn6+uZ1cYDbivqPexn1H/cy6j/uZfb67+7unqvvc1vCfL58+XT9+nX5+PjIx8fndrwlcMdITU1VamqqvL29qf+451D/cS+j/uNeRv3Hvcxe/3Ozi73EAHgAAAAAAFgOYR4AAAAAAIu5LWHey8tLvXv3lpeX1+14O+COQv3HvYz6j3sZ9R/3Muo/7mW3q/7fltHsAQAAAACA69DNHgAAAAAAiyHMAwAAAABgMYR5AAAAAAAshjAPAAAAAIDFEOYBAAAAALAYj1tZadeuXRo/frx27Nih5ORklStXTp07d1ZoaGiWt5GamqrvvvtO33//vY4fPy4fHx9Vr15d/fv3V1BQ0K0UC7gtclr/w8PDtXLlSm3btk2xsbG6cuWKAgMD1bBhQ/Xo0UOFChXK5T0AcsYV/wPSSk5OVuvWrbVv3z6VKVNGv/32m4tLDLiOq+p/fHy8ZsyYod9//13R0dHy9PRUqVKl1LBhQ/Xu3TuXSg/kjCvq/6VLlzRz5kytXLlSMTEx8vLyUsmSJfXiiy+qTZs2yp8/fy7uAXBrfvnlF0VERGj37t06cOCAkpOT9fnnn6tly5bZ2o6rM3C2H023detWdevWTZ6enmratKkKFiyo33//XTExMerfv7969uyZpe18+OGH+uGHH1SuXDk9/fTTOnv2rJYtW6b8+fNr/vz5KleuXLZ3Bshtrqj/tWvX1vnz5/X4448rJCREbm5u2rZtm/bs2aOHHnpI8+fPV7FixW7D3gDZ56r/AWmNHTtWs2bN0pUrVwjzuKO5qv7//fff6ty5s6Kjo1WrVi2FhIQoKSlJJ06c0N9//61ff/01l/cEyD5X1P9Lly6pZcuWio6O1uOPP67HHntMSUlJWr9+vU6cOKEaNWpo5syZypePzsO4szRo0ECxsbEqUqSIfHx8FBsbe0th3uUZ2JYNycnJtkaNGtkeffRR219//WXMv3z5sq1p06a2Rx55xHb06NGbbmfLli228uXL2zp27GhLTEw05m/evNkWHBxse+mll7JTLOC2cFX9nzJliu3UqVOmeampqbaPP/7YVr58edv//vc/VxcdcAlXnQNp7d692/bII4/YZs+ebStfvrzt2WefdXGpAddwVf2/fv26rVWrVrZKlSrZtmzZku77AHcaV9X/qVOn2sqXL28bOnSoaX5iYqKtVatWtvLly9u2bdvm6uIDObZp0yZbTEyMzWa78V2+fPnytp9++ilb28iNDJyty15hYWE6ceKEXnjhBT3yyCPGfD8/P/Xq1UvXr1/XwoULb7qdBQsWSJL69esnLy8vY37NmjVVp04dbd++XUePHs1O0YBc56r63717d5UoUcI0z83NTb169ZIkbd++3bUFB1zEVeeAXVJSkgYOHKjHHntML7/8cm4UGXAZV9X/FStW6M8//1TXrl1Vo0YNp+UeHrd0BySQq1xV/6OjoyVJTz/9tGm+l5eXateuLUk6e/asC0sOuEatWrUUGBiYo23kRgbOVpjftm2bJKlOnTpOy+wnoP01mdm6dat8fHxUtWpVp2X2bRNocKdxVf3PiP0LnLu7+y1vA8hNrj4HJkyYoOPHj+uzzz6Tm5ubawoJ5BJX1f9ly5ZJkp577jmdPHlS8+bN09SpU7V8+XIlJCS4sMSA67iq/v/73/+WJG3YsME0Pzk5WZs3b1aBAgVUpUqVnBYXuCPlRgbO1uXfY8eOSZJKly7ttMzf319FihTR8ePHM93GlStXdPr0aZUvXz7d0GK/8d/+XsCdwhX1PzM//fSTpP/7pwjcaVx5DuzatUvTp09X//79VaZMGVcWE8gVrqr/u3fvliRFRETo888/V1JSkrGsaNGiGjNmjKpXr+6aQgMu4qr636ZNG/3yyy+aMWOGdu/erUcffVTJycnasGGDLl68qJEjR+q+++5zdfGBPJdbGThbLfPx8fGSpIIFC6a73M/PT5cvX850G/blfn5+GW4j7XsBdwpX1P+M7N27VxMnTlSxYsX02muv3XIZgdzkqnMgKSlJ77//vkJCQtS1a1eXlhHILa6q//YuxEOGDFHnzp21bt06bdmyRR988IEuX76sN998U3Fxca4rOOACrqr/BQoU0Jw5c9SsWTNt27ZNM2bM0Jw5c4wu/Om1WAJ3g9zKwAwVCeSx6Oho9ejRQykpKRo1apSKFi2a10UCctWYMWN0/PhxDR06lNtKcM+x/f+HCNWrV08DBgzQ/fffr6JFi6pTp07q0qWLLl++rB9//DGPSwnkjnPnzunVV1/Vzp07NXXqVIWHh2vTpk365JNPtHDhQrVt21YXL17M62IClpGtMG+/YpDRlbf4+PgMr9jZ2ZdndNXBPj+jqxZAXnFF/XcUGxurzp0769y5cxo3bly6gyEBdwpXnAN//fWXZs2apZ49eyo4ONjlZQRyi6v+B9i306BBA6dl9evXl/R/XfGBO4Wr6v+wYcMUFRWlcePG6emnn1bBggUVEBCgtm3b6t1331V0dLS++eYbl5YduBPkVgbOVpi39+VP756Yixcv6vz58+neS5OWj4+PihcvrpiYGKWkpDgtt98nYH8v4E7hivqfVkxMjDp16qS4uDiNGTPG+BIH3KlccQ7s379fKSkpGj9+vIKDg00/knT06FEFBwfriSeecHn5gZxw1f8A+xgRhQoVclpmn5eYmJiDkgKu56r6v27dOhUuXFgPP/yw0zJ7g8Zff/2Vs8ICd6DcysDZCvNPPvmkJGnjxo1OyzZt2iRJqlat2k23U61aNV25ckWRkZFOy+zbtr8XcKdwVf2XbgT5V155RXFxcRo9erQaNWrkuoICucQV50BQUJBat26d7o9048p169at1aJFC9cWHsghV/0PsAeWQ4cOOS2zz8vp448AV3NV/U9KSlJ8fLxp4Ee7c+fOSZLpkV3A3SQ3MnC2wnzNmjVVqlQpLVmyRHv37jXmx8fHa9KkSfLw8NCLL75ozD937pwOHz5snJx2bdu2lXTjvsm0J/OWLVu0ceNGPfnkk4xujDuOq+q/PcifOnVKo0aN0jPPPHPb9gHICVecA1WrVtVnn32W7o8kBQQE6LPPPtMHH3xw+3YMyAJX/Q9o2bKlvLy89O233+rUqVOm7UyZMkWS9Pzzz+fy3gDZ46r6X7VqVV2/fl2TJk0yzU9KSjLm8TQHWN3tzMBuNvtILFkUFham1157TZ6ennrhhRfk5+en33//XTExMerXr5/eeOMN47Xjx4/XhAkT1Lt3b/Xp08e0nQ8++EALFixQuXLl9PTTT+vs2bNatmyZ8ufPr/nz56tcuXLZ2hHgdnBF/W/QoIFiY2NVuXLldJ/XKsnpfAHuFK76H5Ce4OBglSlTRr/99ltu7gJwy1xV/+fMmaMhQ4aocOHCeuaZZ+Tl5aW1a9cqNjZW7dq10+DBg2/3rgE35Yr6v3fvXr300ktKSEhQpUqVVLVqVSUmJmrjxo2Kjo5WhQoVNG/ePOXPnz8vdhHI0IIFCxQRESFJOnDggP766y9VrVrVuL2kUaNGRk/b25mBs/WceelG97DvvvtO48aN0/Lly5WcnKxy5cqpb9++atasWZa3M3jwYAUHB+v777/XnDlz5OPjo/r16/PMYdzRXFH/Y2NjJUk7duzQjh070n0NYR53Klf9DwCsyFX1v1OnTgoMDNTXX3+tpUuXKiUlReXKlVPPnj2NlhvgTuOK+h8SEqKFCxdqypQpCgsL09y5c+Xu7q6HHnpIffr0Ubdu3QjyuCNFRERo0aJFpnmRkZFGl/nAwMAs3Tbr6gyc7ZZ5AAAAAACQt3jOPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACLIcwDAAAAAGAxhHkAAAAAACyGMA8AAAAAgMUQ5gEAAAAAsBjCPAAAAAAAFkOYBwAAAADAYgjzAAAAAABYDGEeAAAAAACL+X9dO7s6NvtI3QAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1000x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot\n",
    "_, ax = plt.subplots(figsize=(10, 4))\n",
    "ax.hist(Ty, np.arange(0, 1.01, 0.05))\n",
    "ax.set(xlim=(0, 1), yticks=[])\n",
    "ax.set_title('Light speed example\\ndistribution of marginal posterior p-values');\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "bda_py_demos",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
